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Separated  Component-Based  Restoration  of 
Speckled  SAR  Images 

Vishal  M.  Patel,  Member,  IEEE,  Glenn  R.  Easley,  Rama  Chellappa,  Eellow,  IEEE  and  Nasser  M. 

Nasrabadi,  Eellow,  IEEE. 


Abstract — Many  coherent  imaging  modalities  such  as  Synthetic 
Aperture  Radar  (SAR)  suffer  from  a  multiplicative  noise  com¬ 
monly  referred  to  as  speckle  which  often  makes  the  interpretation 
of  data  difficult.  An  effective  strategy  for  speckle  reduction  is 
to  use  a  dictionary  that  can  sparsely  represent  the  features  in 
the  speckled  image.  However,  such  approaches  fail  to  capture 
important  salient  features  such  as  texture.  In  this  paper,  we 
present  a  speckle  reduction  algorithm  that  handles  this  issue 
hy  formulating  the  restoration  problem  so  that  the  structure  and 
texture  components  can  be  separately  estimated  with  different 
dictionaries.  To  solve  this  formulation  an  iterative  algorithm 
based  on  surrogate  functionals  is  proposed.  Experiments  indicate 
the  proposed  method  performs  favorably  compared  to  state-of- 
the-art  speckle  reduction  methods. 

Index  Terms — Synthetic  Aperture  Radar,  speckle,  multiplica¬ 
tive  noise,  image  restoration. 

I.  Introduction 

Coherent  imaging  systems  such  as  SAR,  holography,  ultra¬ 
sound,  and  synthetic  aperture  sonar  suffer  from  a  multiplicative 
noise  known  as  speckle  [1].  Speckle  appears  when  objects 
illuminated  by  coherent  radiation  have  surface  features  that 
are  rough  compared  with  the  illuminating  wavelength.  It  is 
caused  by  the  constructive  and  destructive  interference  of 
the  coherent  returns  scattered  by  many  elementary  reflectors 
within  the  resolution  cell.  Speckle  can  make  the  detection 
and  interpretation  difficult  for  automated  as  well  as  human 
observers.  In  some  cases,  it  maybe  important  to  remove 
speckle  to  improve  applications  such  as  compression,  target 
recognition  and  segmentation. 

Many  algorithms  have  been  developed  to  suppress  speckle 
noise  [2],  [3],  [4],  [5],  [6],  [7].  One  of  the  simplest  approaches 
for  speckle  noise  reduction  is  known  as  multi-look  processing. 
It  involves  non-coherently  summing  the  independent  images 
formed  from  L  independent  pieces  of  the  phase  history.  The 
averaging  process  reduces  the  noise  variance  by  a  factor  of 
L.  However,  this  often  results  in  the  reduction  of  the  spatial 
resolution.  Other  types  of  speckle  reduction  methods  are  based 
on  spatial  local  filtering  performed  after  the  formation  of  the 
SAR  image.  Various  filters  have  been  developed  that  avoid  the 
loss  in  spatial  resolution  [2],  [8],  [9],  [10],  [11].  Some  of  these 
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methods  are  based  on  a  window  processing  of  the  noisy  image. 
Consequently,  their  performance  depends  significantly  on  the 
type,  direction  and  the  size  of  the  filter  used.  Furthermore, 
some  of  these  filters  often  fail  to  preserve  sharp  features  such 
as  edges. 

To  overcome  some  of  these  limitations,  wavelet-based  meth¬ 
ods  are  often  utilized  [12],  [13],  [14],  [15],  [16],  [17]  in  which 
noise  shrinkage  is  applied  to  the  detail  wavelet  coefficients  of 
the  noisy  image  [18], [19].  Since  speckle  is  multiplicative  in 
nature,  some  of  these  methods  often  apply  logarithm  transform 
to  SAR  images  to  convert  the  multiplicative  noise  into  additive 
noise.  After  applying  soft  or  hard  thersholding  to  the  wavelet 
coefficients  of  the  logarithmically  transformed  image,  an  ex¬ 
ponential  operation  is  employed  to  convert  the  logarithmically 
transformed  image  back  to  the  original  multiplicative  format. 

It  is  well  known  that  shrinkage-based  denoising  algorithms 
rely  on  the  sparsity  of  the  representation.  A  fixed  transform 
such  as  a  wavelet  transform  can  represent  a  piecewise  smooth 
image  sparsely  but  it  may  also  fail  to  represent  an  image 
with  textures  sparsely.  As  a  result,  the  overall  denoising 
performance  of  a  fixed  transform  on  an  image  containing  both 
piecewise  smooth  and  texture  components  can  be  inadequate. 

Another  popular  approach  for  restoring  speckled  images 
involves  total  variation  (TV)  regularization  [20]  where  the 
underlying  reflectivity  image  is  assumed  to  be  piecewise 
smooth  [21],  [22],  [23].  It  has  been  shown  that  TV  regular¬ 
ization  often  yields  images  with  the  stair-casing  effect  [24]. 
As  a  result  the  estimated  image  usually  contains  constant 
regions  and  fine  details  such  as  textures  are  often  removed. 
To  deal  with  some  of  these  effects,  a  hybrid  method  that 
uses  coefficient  thresholding  and  TV  regularization  on  the 
logarithm  of  the  magnitude  image  or  the  log-image  domain 
was  recently  proposed  in  [25].  In  particular,  hard-thresholding 
on  the  curvelet  coefficients  of  the  log-image  is  first  applied. 
Then  a  variational  method  that  uses  an  fidelity  to  the 
thresholded  coefficients  and  a  TV  regularization  in  the  log- 
image  domain  is  applied. 

Several  methods  have  been  proposed  that  use  a  combined 
dictionary  approach  to  image  restoration.  Suppose  that  we  are 
given  M  different  dictionaries  D^,  rn  =  1,...,M,  then  one 
can  obtain  M  different  estimates  of  x  by  applying  either  hard 
or  soft  threshold  to  the  coefficients  from  each  corresponding 
dictionary.  Let  be  the  resulting  estimate  from  the  ith 
dictionary.  Then,  a  simple  estimator  of  x  is  given  by  averaging 
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M  individual  estimates 


X 


This  is  essentially  the  idea  behind  translation-invariant  wavelet 
denoising  by  cycle  spinning  [26],  where  the  final  estimate 
is  obtained  by  averaging  the  estimates  from  the  thresholded 
orthogonal  wavelet  coefficients  of  translated  versions  of  the 
original  image.  This  method  has  been  applied  for  speckle 
reduction  in  SAR  imagery  [12].  This  simple  method  suffers 
from  some  issues  in  practice  as  it  weights  equally  both 
good  and  bad  quality  estimates.  To  deal  with  this  problem, 
a  Bayesian  framework  to  optimally  combine  the  individual 
estimators  was  proposed  in  [27].  This  method  weights  each 
estimate  at  each  sample  according  to  the  significance  that 
the  elements  in  the  dictionary  have  in  synthesizing  x^  at 
the  same  sample.  This  method  is  effective,  however,  it  can  be 
very  time  consuming.  A  similar  approach  was  also  proposed 
in  [28]. 

A  novel  cartoon  and  texture  image  component-based 
restoration  for  SAR  images  was  proposed  in  [29].  Specifically, 
a  SAR  image,  considered  as  a  function  /,  is  to  be  decomposed 
into  a  sum  of  two  components  f  =  u  +  v,  where  u  represents 
the  cartoon  or  geometric  (i.e.  piecewise  smooth)  components 
of  /,  and  V  represents  the  oscillatory  or  textured  components. 
The  second  component  essentially  accounts  for  noise  and  the 
texture  elements.  In  the  proposed  technique,  u  is  estimated 
and  is  considered  as  the  restored  (despeckled)  image  while  the 
textured  components  of  v  are  not  attempted  to  be  recovered. 
This  is  problematic  since  discarding  the  texture  components 
may  result  in  the  loss  of  important  salient  features  in  a  SAR 
image.  Figure  1  provides  an  example  of  how  an  image  sep¬ 
arates  into  the  structural  and  textured  components  indicating 
the  importance  of  retaining  these  textural  elements. 

Motivated  by  recent  advances  in  sparse  representation- 
based  image  separation  [30],  we  propose  a  similar  separation- 
based  despeckling  method  so  that  the  image  is  decomposed 
into  a  sum  of  piecewise  smooth  and  textured  elements.  Our 
formulation  is  based  on  finding  sparse  representations  of  these 
elements  from  dictionaries  specifically  suited  to  compress 
them  and  differs  significantly  by  not  just  treating  the  texture 
and  noise  components  as  complements  of  the  cartoon-based 
estimated  image.  By  taking  advantage  of  the  ability  of  sparse 
representations  in  our  scheme  to  estimate,  we  are  able  to  retain 
important  salient  and  textured  features  in  the  final  estimated 
image. 


A.  Organization  of  the  paper 

Section  II  reviews  the  statistical  models  for  SAR  images. 
In  Section  III,  we  present  the  component  separation  algorithm 
based  on  surrogate  functionals.  In  Section  IV,  we  show  some 
of  the  results  on  both  real  and  simulated  data,  and  present  the 
concluding  remarks  in  Section  V. 

II.  Statistical  Models  for  SAR  Images 

Let  Y  G  be  the  observed  image  intensity,  X  G 

j-jjg  j^Qjgg  fj-gg  image,  and  F  G  be  the  speckle 


Fig.  1.  Image  separation,  (a)  Original  image,  (b)  Structural  components, 
(c)  Textural  components,  (d)  Structural  +  Textural  components.  Note  that  the 
textural  components  in  (c)  contains  some  important  features. 


noise.  Assuming  the  SAR  image  represents  an  average  of  L 
looks,  Y  is  related  to  X  by  the  following  multiplicative  model 

[31] 

y  =  FX,  (1) 

where  F  is  the  normalized  fading  speckle  noise  random 
variable.  It  follows  a  Gamma  distribution  with  unit  mean  and 
variance  and  has  the  following  pdf 

F  >  0,  L  >  1  where  r(.)  denotes  the  gamma  function. 

The  natural  logarithmic  transformation  converts  the  multi¬ 
plicative  noise  model  (1)  into  an  additive  model 

Y  =  ln(y)  =  ln(F)  +  ln(X) 

=  FfX,  (3) 

where  F  =  ln(F)  and  X  =  ln(X).  The  pdf  of  the  random 
variable  F  is  given  by  [32],  [31] 

The  mean  of  F  is  given  by 

F[F]=t^(0,L)-ln(L) 
and  variance  is  given  by 

var(F)  = 

where 

=  inr(f) 
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is  known  as  the  fcth  polygamma  function. 

Note  that  there  exist  many  statistical  models  for  SAR 
images  [33],  some  based  on  mathematical  approaches  and 
others  based  on  physical  approaches.  In  this  paper,  we  use 
the  logarithmic  transform  to  convert  the  multiplicative  noise 
to  an  additive  noise  and  take  into  account  the  appropriate 
distribution  model  to  handle  the  noise. 

III.  Image  Component  Estimation 

Let  y,  f ,  and  x  be  lexicographically  ordered  vectors  of  size 
iV^  representing  Y,F,  and  X,  respectively.  We  assume  that 
the  SAR  reflectivity  field  is  a  superposition  of  two  different 
signals 

x  =  Xp  +  xt,  (5) 

where  Xp  is  the  piecewise  smooth  component  of  x  and  Xf  is 
the  texture  component  of  x.  So  the  additive  model  (3)  can  be 
written  as 

y  =  x  +  f 

=  Xp  +  xt  +  f .  (6) 

We  further  assume  that  Xp  is  compressible  in  a  dictionary 
represented  in  a  matrix  form  as  Dp,  and  similarly,  Xf  is 
compressible  in  a  dictionary  represented  in  a  matrix  form  as 
Dt.  Given  M^,Mt  >  N"^,  the  dictionary  Dp  G 
and  Dj  G  are  chosen  such  that  they  provide  sparse 

representations  of  piecewise  smooth  and  texture  contents, 
respectively.  That  is,  we  assume  there  are  coefficient  vectors 
ttp  G  and  at  G  so  that  Xp  =  DpCip  and 

Xt  =  Dtat.  The  compressibility  assumption  means  that  when 
the  coefficients  are  ordered  in  magnitude,  they  decay  rapidly. 
The  texture  dictionary  Dt  needs  to  contain  atoms  that  are 
oscillatory  in  nature  such  as  those  found  in  the  discrete 
cosine/sine  transform  and  the  Gabor  transform.  The  dictionary 
Dp  should  be  able  to  process  images  with  geometric  features 
such  as  edges.  The  matrix  Dp  could  represent  some  type  of 
wavelet,  shearlet,  curvelet,  or  contourlet  dictionary. 

One  can  recover  the  SAR  reflectivity  field  x  by  estimating 
the  components  Xp  and  xt  via  ap  and  at  by  solving  the 
following  variational  problem: 

ap,at  =  arg  min  A||Q;p||i  +  A||Q;t||i  +^TV{'Dpap) 

Ctp,0't 

+  ^lly-Dpap-Dtat||2,  (7) 

where  TV  is  the  total  variation  (i.e.  sum  of  the  absolute 
variations  in  the  image)  and  for  an  A^-dimensional  vector  x, 
|j.||q  denotes  the  f^-norm,  0  <  q  <  oo,  defined  as 

/  A 

\i=l 

The  two  components  are  the  corresponding  representations 
of  the  two  parts  and  can  be  obtained  by  Xp  =  Dpdp  and 
Xt  =  Dtdt.  This  notion  of  separating  a  signal  into  different 
morphologies  using  sparse  representations  is  often  known  as 
Morphological  Component  Analysis  (MCA)  [30]. 

Instead  of  seeking  the  sparse  sets  of  coefficients  ap  and  at 
directly  and  then  inverting  the  representations,  it  is  possible 


to  directly  seek  the  images  whose  transform  coefficients  or 
dictionary  representations  are  sparse.  This  corresponds  to 
solving  the  following  optimization  problem: 

Xp,Xt  =  arg  min  A||D],Xp||i  +  A||Djxt||i  +-fTV{xp) 

Xp  ,Xt 

+  ^lly-xp-xt||2,  (8) 

where  D|  denotes  the  Moore-Penrose  pseudo  inverse  of  D*. 
Here,  we  have  assumed  that  the  two  redundant  dictionaries  are 
of  full  rank  and  we  can  obtain  the  analysis  operator  from  the 
synthesis  by  using  the  following  relations 

ap  —  Dj^Xp 

at  =  D|xt. 

One  of  the  major  advantages  of  using  (8)  is  that  it  requires 
searching  lower  dimensional  vectors  rather  than  longer  di¬ 
mensional  representation  coefficient  vectors.  This  increases 
numerical  efficiency  and  decreases  memory  constraints. 


A.  Iterative  Shrinkage  Algorithm 

Various  methods  can  be  used  to  obtain  the  solution  of  (7) 

[34]  [35].  In  this  section,  we  derive  a  fast  convergent  iterative 
shrinkage  algorithm  by  a  method  of  using  Separable  Surrogate 
Functionals  (SSF)  to  solve  the  separation  problem  posed  in  (7) 

[35] ,  [36],  [37].  For  simplicity,  we  assume  that  D  =  [Dp,  D*] 
and  discard  the  TV  component  for  the  discussion  given  here. 
The  objective  function  in  (7)  can  then  be  re-written  as 

/(a)  =A||a||i-f  i||y-Da||i  (9) 

where  a  contains  both  the  piecewise  smooth  and  texture  parts. 
Let  ^ 

d(a,ao)  =  |||a  -  ao||2  “  “  ^ao\\l,  (10) 

where  ag  is  an  arbitrary  vector  of  length  N'^  and  the  parameter 
c  is  chosen  such  that  d  is  strictly  convex.  This  constraint  is 
satisfied  by  choosing 

c>  ||D^D||2  =A^ax(D^D), 


where  Amax(D^D)  is  the  maximal  eigenvalue  of  the  matrix 

D^D. 

Adding  (10)  to  (9)  gives  the  following  surrogate  function 

/(a)  =  X\\a\\i  +  hy-T>a\\l  +  ^\\a-ao\\l-h'Da-I)ao\\l. 

(11) 

This  surrogate  function  f{a)  can  be  re-expressed  as 

/(a)  =  A-f-||a||i-f  ijja-xof,  (12) 

c  2 

where  ^ 

xo  =  -D^(y  -  Dao)  +  cto 
c 


and  A  is  some  constant.  Let  (a)+  denote  the  function 
max(a,  0)  and  sign(x)  be  the  signum  function  defined  as 


[  -1  ifa;<0 

signjcc)  =  \  0  for  a;  =  0 
1  1  for  a;  >  0. 
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Given  that 

5a(x)  =  sign(x)(|x|  -  A)+  (13) 

is  the  element-wise  soft-thresholding  operator  with  threshold 
A,  the  global  minimizer  of  the  surrogate  function  is  given  by 


^sol  —  ^X/c  (^o) 

=  ‘5a/c  (^^D’^(y-Dao) +ao^  •  (14) 

It  was  shown  in  [36]  that  the  iterations 

«'=+'=  5A/e(^iD^(y-Da'=)+a'=^  (15) 

converge  to  the  minimizer  of  the  function  /  in  (9)  where  the 
superscript  k  indicates  that  it  is  the  value  for  the  fcth  iteration. 
By  breaking  the  above  iteration  into  the  two  representation 
parts  and  considering  the  TV  term,  we  get  the  following 
iterative  updates  that  essentially  solves  (7) 

=  5A/e  (y  -  (16) 

a^+i  =  (H^Dp5*+i)  (17) 

=  5a/c  (y  -  ,  (18) 

where  H  is  the  undecimated  Haar  wavelet  dictionary.  A  de¬ 
tailed  description  of  the  undecimated  Haar  wavelet  transform 
can  be  found  in  [38].  We  have  replaced  the  TV  correction 
term  by  a  redundant  Haar  wavelet-based  shrinkage  estimate  as 
this  seems  to  give  the  best  results.  This  adjustment  is  applied 
only  to  the  piecewise  smooth  component  to  control  the  ringing 
artifacts  near  the  edges  caused  by  the  oscillations  of  the  atoms 
in  the  dictionary  Dp.  The  same  adjustment  was  used  in  [30] 
and  the  substitution  was  partially  motivated  by  observing  the 
connection  between  TV  and  the  Haar  wavelet  given  in  [39]. 

The  iterations  presented  above  can  be  extended  to  handle 
the  analysis  formulation  in  equation  (8).  This  is  simply  done 
by  modifying  iterations  (16),  (17)  and  (18)  as  follows 

-  Dp.5A/e  (^Dj(y  -  x^  -  xf)  -f  Dtx^^  (19) 
x*+i  =  (H^x^+i)  (20) 

xf+i  -  D*.5a/c  (y  -  x^  -  xf )  +  Djxf)  .  (21) 

We  summarize  the  algorithm  for  recovering  the  two  separated 
components  of  a  SAR  image  in  Fig.  2.  In  step  3  of  the  algo¬ 
rithm  in  Fig.  2,  ||.||oo  denotes  the  foo-norm.  For  an  N  dimen¬ 
sional  vector  X,  it  is  defined  as  |jx||oo  =  max(|a;i|,  •  •  •  ,  |a;Ar|). 

Once  the  two  denoised  components  of  x  are  estimated,  we 
obtain  the  final  estimate  of  x  as 

X  =  exp(xp -f  xt).  (22) 


Input:  y,  c. 

Initialization:  Initialize  fc  =  1  and  set 

x[!  =  0,  X?  =  0,  r°  =  y  —  x°  —  X?,  and  A°  = 

HlIDjylloo  +  IID^yll^). 

repeat: 

1.  Update  the  estimate  of  Xp  and  xt  as 

=  Dp.^A. 

x^  =  (H^ip) 

2.  Update  the  residual  as 


3.  Update  the  shrinkage  parameter  as 

A'=  =  i(||Djr'=||o.  +  ||Dfr'=||o.). 

until:  stopping  criterion  is  satisfied. 

Output:  The  two  components  Xp  =  Xp  and  Xt  =  x*. 


Fig.  2.  The  SSF  iterative  shrinkage  algorithm  to  solve  (8). 


IV.  Experimental  Results 

In  this  section,  we  present  the  results  of  our  proposed 
despeckling  algorithm  and  compare  them  with  the  enhanced 
Lee  filter  [2]  and  some  recent  state-of-the-art  methods  [25], 
[21],  [22],  [30].  We  also  compare  our  results  with  a  Stein- 
Block  thresholding  (SBT)  method  proposed  in  [40].  This 
method  was  shown  to  be  nearly  minimax  over  a  large  class 
of  images  in  the  presence  of  additive  bounded  noise.  This 
method  requires  a  threshold  parameter  which  we  set  to  the 
theoretical  value  4.505  as  derived  in  [40].  Furthermore,  we 
compare  the  performance  of  our  combined  dictionary-based 
approach  to  despecking  with  that  of  a  fixed  transform-based 
despecking  method.  In  particular,  we  apply  soft-thresholding 
on  the  subband  coefficients  of  the  wavelet  transform.  We  call 
the  resulting  method  wavelet-based  thresholding  (WT).  For  the 
MCA  method  [30],  we  use  the  curvelet  transform  to  represent 
the  piecewise  smooth  component  and  2D-DCT  to  representant 
the  texture  component. 

In  Fig.  3,  we  display  the  test  images  used  for  different 
experiments  in  this  paper.  In  these  experiments,  we  use  the 
relative  error  (RE)  and  the  equivalent  number  of  looks  (ENL) 
to  measure  the  performance  of  the  routines  tested: 


ENL  = 


mean^ 
variance  ’ 


Since,  the  logarithmic  transformation  introduces  a  bias  on 
the  final  estimate,  we  correct  it  along  with  the  exponential 
transformation  [12],  [25]. 


where  the  mean  and  the  variance  are  measured  within  a 
homogeneous  region  (see  [21]  and  [2]).  A  large  ENL  value 
corresponds  to  better  speckle  reduction. 
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(C) 


Fig.  3.  Images  used  in  this  paper  for  different  experiments,  (a)  256  X  256 
Cameraman  image,  (b)  512  X  512  Fields  image,  (c)  512  X  512  Nimes  image. 


A.  Dictionaries 

In  our  experiments,  we  use  a  dictionary  corresponding  to 
the  shearlet  transform  to  represent  the  piecewise  smooth  com¬ 
ponent.  Shearlets  extend  the  traditional  wavelets  by  allowing 
for  waveforms  to  be  defined  not  only  at  various  scales  and 
locations  but  also  at  various  orientations  [41],  The  shearlet 
transform  is  best  suited  for  representing  images  with  edges  and 
anisotropic  structures.  Moreover,  numerical  results  give  evi¬ 
dence  to  the  superior  behavior  of  shearlet-based  decomposition 
algorithms  when  compared  to  curvelet-based  and  contourlet- 
based  algorithms  [41],  [42].  This  is  the  reason  why  we  choose 
shearlets  instead  of  curvelets  and  contourlets  in  our  approach. 
A  brief  discussion  about  the  shearlet  transform  is  given  in  the 
appendix.  Fig.  4,  shows  some  atoms  from  a  shearlet  dictionary. 
In  our  implementation,  we  used  the  nonsubsampled  shearlet 
transform  with  the  decomposition  structure  [3,  3,4,4]  which 
determines  the  number  of  directions  in  the  scales  from  coarse 
to  fine.  As  a  result,  the  size  of  Dp  is  Af2  X  56Ar2. 

The  discrete  cosine  transform  is  known  to  closely  approx¬ 
imate  the  Karhunen-Loeve  transform  for  a  class  of  random 
signals  known  as  first-order  Markov  processes  that  model 
several  real-world  images.  Explicitly  the  coefficients  of  the 
DCT  for  w  N  X  N  image  x  =  {xm,n)  are 

2  N-lN-l 

Xj,k  =  X!  X!  cos(7r(m  -  l/2){j  -  l/2)/N) 

m—0  n—0 

X  cos(7r(n  —  l/2)(fc  —  l/2)/N). 

It  is  clear  the  oscillatory  nature  of  the  basis  functions  can 
represent  repetitious  elements  common  to  texture  components 
which  is  why  they  are  commonly  use  for  this  purpose.  A  few 
atoms  from  the  2D-DCT  dictionary  are  shown  in  Fig.  5.  The 
number  of  atoms  in  the  dictionary  Dt  are  64  x  N"^. 


Fig.  4.  A  few  atoms  from  a  shearlet  dictionary.  Each  block  represents 
the  result  of  the  applying  the  shearlet  transform  for  a  particular  scale  and 
orientation  after  applying  it  to  a  centered  impulse  response. 


Fig.  5.  A  few  atoms  from  the  2D-DCT  dictionary.  Each  block  represents  the 
result  of  applying  the  inverse  2D-DCT  transform  to  impulse  responses  which 
are  centered  at  various  locations. 


B.  Parameters 

From  the  discussion  in  section  III,  the  parameter  c  should 
be  chosen  such  that  c  >  Amax(DpDj  4-DtD^).  This  can  be 
satisfied  by  choosing  c  >  2.  In  particular,  the  value  we  choose 
was  c  =  3. 

The  Haar  shrinkage  value  7^  in  step  1  of  the  algorithm 
is  3(tJ  ,  where  is  the  standard  deviation  of  the  noise 

xLp  ‘  iLp 

estimated  by  using  a  median  estimator  on  the  finest  scale  of 
the  Haar  wavelet  coefficients  of  Xp  [19]. 


C.  Stopping  rule 

We  change  the  threshold  value  of  during  each  iteration 
according  to 

A'=  =  i(||Djr'=||o.  +  ||Dfr'=||o.) 

and  stop  the  iterations  when  A^  <  Ttr,  where  T  w  2.1  [43]. 
This  way  we  stop  the  iterations  when  the  residual  is  at  the 
noise  level  and  the  noise  is  rejected  in  each  component.  Note 
that  this  step  of  the  algorithm  assumes  the  noise  variance  to 
be  known.  This  is  not  a  problem  since  the  noise  variance  can 
easily  be  estimated  by  using  a  median  estimator  on  the  finest 
scale  of  the  wavelet  coefficients  of  y  [19]. 


D.  Results  on  simulated  data 

In  Table  I,  we  report  the  results  of  experiments  on  the 
simulated  data  with  various  number  of  looks.  We  did  not 
have  access  to  the  codes  of  [25],  [21]  and  [22].  Hence,  we 
report  the  relative  error  values  reported  in  the  corresponding 
papers.  Figures  6  and  7  show  the  noisy  and  restored  images 
for  some  of  the  experiments  with  the  simulated  data.  It  can 
be  seen  from  these  figures  and  the  results  in  Table  I  that 
our  method  performs  favorably  over  some  of  the  competitive 
methods  for  speckle  reduction.  In  particular,  this  is  the  case 
when  the  number  of  looks  is  greater  than  three.  Furthermore, 
these  results  clearly  indicate  that  an  improvement  is  achieved 
when  a  combined  dictionary  approach  is  used  to  restore  a  SAR 
image  as  can  be  seen  by  comparing  the  results  of  our  method 
with  that  of  SET,  MCA  and  WT  in  Table  I. 


(c)  (d) 


Fig.  6.  (a)  Noisy  image,  L  =  4,  RE  =  0.498.  (b)  Restored  image  using  our 
method,  RE  =  0.118.  (c)  Noisy  image,  L  =  A,  RE  =  0.500.  (d)  Restored 
image  using  our  method,  RE  =  0.065. 

Using  the  Cameraman  image,  in  Figure  8  and  Figure  9, 
we  show  the  evolution  of  the  objective  function  and  relative 
error,  respectively,  as  we  vary  the  number  of  looks.  Note  that 
in  Figure  9,  the  relative  error  decreases  significantly  after  five 
iterations  and  saturates  around  the  tenth  iteration,  showing  that 
the  proposed  method  is  efficient  and  requires  less  number  of 
iterations  compared  to  [21]. 

E.  Results  on  real  SAR  data 

In  the  second  set  of  experiments,  we  use  the  real  SAR 
images  shown  in  Figure  10(a)  and  Figure  10(b).  These  images 
were  collected  using  the  Sandia  National  Laboratories  Twin 
Otter  SAR  sensor  payload  operating  at  X  band.  Since  the  true 
reflectivity  fields  are  not  available,  we  use  ENL  to  measure  the 


(e)  (f) 

Fig.  7.  Experiment  with  the  Nimes  image,  (a)  Noisy  image,  L  =  10,  RE  = 
0.315.  (b)  Restored  image  using  our  method,  RE  =  0.163.  (c)  Noisy  image, 
L  =  A,  RE  =  0.501.  (d)  Restored  image  using  our  method,  RE  =  0.207. 
(e)  Noisy  image,  L  =  1,  RE  =  1.001.  (f)  Restored  image  using  our  method, 
RE  =  0.290. 


performance  of  our  method.  The  value  of  ENL  is  estimated 
from  the  two  32  x  32  homogeneous  regions  (shown  with  white 
boxes).  We  refer  to  the  region  left  side  of  the  image  as  R1  and 
the  right  side  of  the  image  as  R2.  The  estimated  ENL  values 
are  reported  in  Table  IT  As  can  be  seen  from  this  table,  our 
method  outperforms  the  WT,  MCA  and  SET  methods. 

The  despeckling  results  from  various  methods  are  shown 
in  Eigures  ll(b)-(d)  and  Eigures  12(b)-(d).  The  ratios  of  the 
original  image  to  the  filtered  images,  referred  to  as  the  noise 
images,  are  shown  in  Eig.  ll(e)-(h)  and  Eigures  12(e)-(h). 
It  is  evident  from  these  figures  the  single  dictionary-based 
reconstructions  such  as  SET  and  WT  suffer  from  noticeable 
artifacts.  The  MCA  method  based  on  curvelets  and  2D- 
DCT  provides  good  reconstruction,  however,  it  removes  a  lot 
of  point  targets.  Our  combined  dictionary  approach  clearly 
provides  good  reconstructions  and  removes  most,  if  not  all,  of 
these  artifacts  and  preserves  point  targets.  Note  that  we  have 
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TABLE  I 

Relative  errors  for  various  experiments. 


Image 

L 

Noisy 

Ours 

SET  [40] 

MIDAL  [21] 

[25] 

[22] 

WT 

Lee 

MCA 

Cameraman 

13 

0.277 

0.083 

0.104 

0.090 

- 

0.098 

0.098 

0.170 

0.090 

Cameraman 

10 

0.315 

0.090 

0.111 

0.097 

0.091 

- 

0.105 

0.171 

0.101 

Cameraman 

4 

0.498 

0.118 

0.144 

0.124 

0.131 

- 

0.135 

0.178 

0.140 

Cameraman 

3 

0.573 

0.129 

0.156 

0.130 

- 

0.151 

0.147 

0.182 

0.151 

Cameraman 

1 

0.990 

0.184 

0.220 

0.167 

0.192 

- 

0.228 

0.211 

0.210 

Fields 

10 

0.316 

0.054 

0.058 

0.056 

0.055 

- 

0.063 

0.063 

0.052 

Fields 

4 

0.500 

0.065 

0.073 

0.066 

0.066 

- 

0.077 

0.804 

0.071 

Fields 

1 

1.000 

0.099 

0.114 

0.089 

0.096 

- 

0.159 

0.135 

0.110 

Nimes 

10 

0.315 

0.163 

0.312 

0.170 

0.174 

- 

0.195 

0.273 

0.186 

Nimes 

4 

0.501 

0.207 

0.223 

0.217 

0.217 

- 

0.247 

0.277 

0.221 

Nimes 

1 

1.001 

0.290 

0.312 

0.301 

0.314 

- 

0.346 

0.298 

0.320 

Fig.  8.  The  objective  function  value  as  a  function  of  iteration  number  for 
the  experiments  with  a  Cameraman  image  with  various  number  of  looks. 


Fig.  9.  Relative  eiTor  vs.  number  of  iterations  curves  for  the  experiments 
with  a  Cameraman  image  with  various  number  of  looks.. 


Fig.  10.  Real  SAR  images  used  for  the  experiments. 

not  compared  our  method  on  the  real  SAR  images  with  the 
other  methods  since  their  implementations  were  not  available 
to  us. 

F.  Computational  efficiency 

In  our  image  separation-based  despeckling  method,  the  most 
computationally  intensive  part  is  in  finding  the  coefficients 
from  a  shearlet  dictionary.  Using  MatLab  on  a  Windows 
system  with  Intel  Core  2  CPU  2.16  GHz/3.00  GB  processor, 
one  iteration  of  the  entire  algorithm  takes  around  4.15  seconds. 
On  average  our  method  takes  about  45  seconds  to  process  an 
image  of  size  256  x  256.  The  performance  of  our  algorithm 
can  be  enhanced  by  using  a  more  efficient  shearlet  transform 
implementation  which  is  parallelizable.  Each  iteration  of  our 
algorithm  has  the  complexity  of  0{N^  log2(iV)). 

G.  Discussion 

Note  that  we  have  applied  our  component  separation  method 
on  the  log-transformed  images.  However,  one  can  directly 
apply  this  method  on  the  SAR  image  without  applying  the  log 
transformation  to  separate  the  piecewise  smooth  component 
and  the  texture  component.  For  instance,  SAR  model  (1)  can 
be  re-written  as 

y  =  xf 

=  X  -h  (f  -  l)x,  (23) 

where  1  =  [1,  •  •  •  ,1]^-  In  this  setting,  the  first  and  the 
second  term  in  (23)  can  be  viewed  as  Xp  and  xj,  respectively. 
In  particular,  if  one  has  designed  dictionaries  specifically  to 


(e)  (f)  (g)  (h) 


Fig.  11.  Despeckled  images.  Restored  using  (a)  the  SET  method,  (b)  the  MCA  method,  (c)  wavelet-based  thresholding  and  (d)  our  method,  (e)-(h)  are  the 
noise  images  corresponding  to  the  filtered  images  in  (a)-(d),  respectively. 


(e)  (f)  (g)  (h) 


Fig.  12.  Despeckled  images.  Restored  using  (a)  the  SET  method,  (b)  the  MCA  method,  (c)  wavelet-based  thresholding  and  (d)  our  method,  (e)-(h)  are  the 
noise  images  corresponding  to  the  filtered  images  in  (a)-(d),  respectively. 
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TABLE  II 

The  ESTIMATED  ENL  VALUES. 


Image 

Region 

Original 

WT 

SET 

Ours 

Lee 

MCA 

Fig.  11(a) 

R1 

3.291 

30.366 

49.028 

79.357 

5.518 

50.23 

Fig.  11(a) 

R2 

3.236 

30.428 

84.936 

133.962 

8.263 

90.89 

Fig.  12(b) 

R1 

4.117 

44.717 

75.041 

146.379 

11.35 

85.75 

Fig.  12(b) 

R2 

3.997 

37.253 

74.641 

101.491 

11.30 

86.70 

compress  these  terms  then,  one  can  view  the  first  term  as  the 
restored  image.  However,  designing  dictionaries  to  compress 
the  individual  terms  in  (23)  is  nontrivial.  One  can  adapt 
dictionary  learning  methods  [44]  to  learn  these  components. 
This  in  turn  would  require  a  collection  of  training  images 
containing  only  the  piecewise  smooth  reflectivity  fields  and 
images  containing  only  signal-dependent  noise.  This  formu¬ 
lation  is  also  relevant  in  the  case  when  the  SAR  images 
contain  a  strong  additive  noise  component  where  the  purely 
multiplicative  model  may  not  be  adequate. 

V.  Conclusion 


for  a  €  K."*",  s  G  R,  and  t  €  The  operator  SH  defined  by 

SHf{a,s,t)  =  (fj-ipast) 


is  referred  to  as  the  continuous  shearlet  transform  of  /  G 
L^(]R).  It  is  dependent  on  the  scale  variable  a,  the  shear  s, 
and  the  location  t. 

The  collection  of  discrete  shearlets  is  given  by 
{fjAk  =  I  det  -  k)  :  j,e  G  Z,k  G 


where 


B  = 


1  1 
0  1 


A  = 


2  0 
0  V2 


We  have  proposed  a  new  method  of  speckle  reduction 
in  SAR  imagery  based  on  separating  an  image  into  various 
components.  Unique  to  this  approach  is  the  ability  to  use 
specific  dictionaries  of  representations  suited  for  separation 
with  an  iterative  scheme  that  is  able  to  retain  important  fea¬ 
tures.  The  experiments  show  this  method  performs  favorably 
compared  to  other  competitive  methods.  This  new  process  is 
also  valuable  for  many  SAR  image  understanding  tasks  such  as 
road  detection,  railway  detection,  ship  wake  detection,  texture 
segmentation  for  agricultural  scenes  and  coastline  detection. 
In  addition,  specific  dictionaries  could  be  designed  to  be  used 
with  this  procedure  to  capture  unique  signatures  while  dealing 
with  the  speckle  removal. 


Appendix 

The  Shearlet  Transform: 

The  shearlet  construction  can  be  considered  as  a  natural 
extension  of  wavelets  into  two-dimensions  [45].  Its  representa¬ 
tive  elements  are  defined  by  the  two-dimensional  affine  system 

{fastix)  =  I  deiMas\~^  -t)  :  tG  R^}, 


where 


is  a  product  of  a  shearing  and  anisotropic  dilation  matrix  for 
(a,  s)  G  X  R.  The  generating  function  ip  is  such  that 


^(C)  =  ^(Ci,6) 


where  ipi  is  a  continuous  wavelet  for  which  iki  G  (^“(R) 
with  supp  C  [—2, 1/2]  U  [1/2,  2],  and  1/2  is  chosen  so  that 
4*2  G  (^“(R),  suppT'2  C  [—1, 1],  with  T'2  >  0  on  (-1,1),  and 
1 1 1/2 II 2  =  1-  Under  these  assumptions,  a  function  /  G  L^(R^) 
can  be  represented  as 

p  pOO  pOO  7 

f{x)=  /  /  /  {f,'>Past)'ipast{x)-^  dsdt, 

J —00  Jo  ^ 


Shearlets  form  a  Parseval  frame  (tight  frame  with  bounds  equal 
to  1)  for  L^(R^)  given  the  appropriate  choice  of  the  generating 
function  ip  (see  [41]  for  details).  An  M-channel  filterbank 
implementation  can  be  done  by  using  the  techniques  given  in 
[46].  As  a  consequence,  its  implementation  has  a  complexity 
of  0{N'^  log2(W))  for  an  N  X  N  image. 
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